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0. INTRODUCTION 

This final report summarizes most of the research sponsored by the 
National Aeronautics and Space Administration under Grant NGR - 33-006- 
020 during the period January 1, 1972 to January 1, 1973. The research 
supported by this grant deals mostly with problems of digital data trans- 
mission and includes some new work on computer communication networks. 
There are four self-contained sections labeled I through IV. Each section 
has its own figures, references and equation numbering. . 

Section I, Signal Processing with Finite State Machines, continues the 
work on finite memory detection. New results for a time -invariant 
machine are given, optimum time -varying machines for detection are pre- 
sented, and the structure and performance of these machines are developed. 
Problems of practical significance are discussed. Many talks have been 
given on this subject. One paper has been accepted for publication, and 
a second is almost completed. Invited papers were presented at the 
IEEE Convention, in New York in March 1972 and at the Southeastern 
Systems Symposium in North Carolina in March 1973. These papers have 
been published in the records of the symposia. Two papers were presented 
at the International Symposium on Information Theory in California in 
January 1972. This work was performed by R. R. Boorstyn, P. F. Lynn, 
and R. W. Muise and is continuing. 

Section II discusses work on Signal Parameter Estimation from Dis- 
crete-Time Observations. New results relevant to the problem of esti- 
mating several single-frequency tones from a finite number of noisy, 
discrete -time observations are presented. This problem has application 
to data systems testing, radar, and other measurements. A paper will be 
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presented at the International Symposium on Information Theory m Israel 
in June 1973 and will be submitted for publication. This work was performed 
by R, R. Boor styn and D. C. Rife. 

Section III contains work on Digital Filtering for Radar Signal Proces- 
sing Applications. A novel approach in synthesizing digital filters for signal 
processing applications is presented. This synthesis method takes advan- 
tage of the known signal waveform structure and results in many fewer 
digital computations as compared to convolution processing. This approach 
is particularly suited to synthesis of matched filters for radar signal proces- 
sing and yields matched or approximately matched filters which simulta- 
neously have very low storage and very low computational requirements. 

A paper has been published in the Transactions on Audio and Electronics, 
IEEE, March 1972. This work was performed by R. R, Boor styn and J. D. 
Echard. 

Section IV contains work On Multiple Server Queues Where Not All 
Servers are Identical. This is an attempt to derive some properties of net- 
works of queues by considering the individual outputs of multiple server 
queues. It is shown under what conditions the outputs retain the Poisson 
character of the input. This work was performed by R. R. Boor styn and 
P. McGregor. 

The PIB faculty participating in this program were R. R. Bo or styn, 
who prepared the final report, R. A. Haddad, M. Schwartz, and J. K. Wolf , 
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I. SIGNAL PROCESSING WITH FINITE STATE MACHINES 

This paper is concerned with the implementation of communications 
receivers using digital processors or digital computers and the effect of this 
implementation on performance. One can consider conventional (analog) re- 
ceivers to consist of filters, summers, multipliers, comparators, matched 
filters, phase-locked loops, etc. The usual method in converting these re- 
ceivers into digital devices is to replace each analog component with its 
digital counterpart. Thus, after the required sampling, digital filters re- 
place the original analog filters, summers and multipliers are replaced by 
digital summers and multipliers, etc. If the sampling is fast enough there 
should be little deleterious effect on performance. 

However, there is one further constraint that is the dominant concern 
of this paper. Every coefficient, every data sample, and the results of 
every operation must be represented and stored by a finite length word, 
i. e. , a finite number of bits. For a standard computer approximately 30 
bits are available for each word and this should be sufficient For special 
purpose computers, or in time- sharing operation where the device is to be 
used simultaneously for many similar operations, it is important to in- 
vestigate the effect of reducing the storage capability considerably. It 
is clear that as the number of bits is reduced, the performance will ulti- 
mately deteriorate. If the amount of storage is increased from this point, 
the performance will be reasonably satisfactory. However, if the size of 
storage is reduced further the question arises; is there an entirely new 
digital structure that will perform well with a minimum amount of storage? 
This question is answered affirmatively in this paper. 

To illustrate, consider a simple binary receiver. One of two equally 
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likely signals are to be sent, and n independent samples X^, X^, • * • » 
are received. If one of the signals (denoted by hypothesis H q ) is sent then 
each of the X, has density f (x). Otherwise (H^ they have density f } (x) . 
The optimum receiver forms 


f ,(X ) 

(1) \(k) = i n y (x ' T + ^ k_1 ^ k = lf 


n, 


where \{o) = 0, and decides in favor of Hj if X(n)>0. This is just a re- 
cursive version of the likelihood ratio. 

In general all the n data samples must be remembered with infinite 
precision to make the optimum decision. If the data is assumed to be col- 
lected sequentially in time, then equation (1) states that after k samples are 
received these first k samples of data need effectively only be "remembered” 
by storing the single number X(k) . Then X(n) is used to make the decision, 
But the X(k) must still be stored with infinite precision. It can be shown 
that, in simple examples, merely uniformly quantizing each X(k) to a finite 
number of bits and using equation(l) results in a deterioration of perfor- 
mance if the number of samples is roughly greater than the number of le- 


vels of quantization. 

There are, however, receiver structures which incorporate funda- 
mentally this finite storage constraint, here called "finite memory", and 


perform almost as well as the optimum infinite memory receiver 


/> 1 r 


equation (1). These receivers can be modelled as finite state machines. 
After the (k- l) St sample the machine may be in one of m states, e. g. , the 
m levels of a log m bit quantizer. This state at time k-1, S^_|> an< ^ ^he 

Li 

next sample X k are used to determine the state at time k. The final de- 
cision is based upon the state at time n. An example of such a finite 
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machine suitable for detection is shown in Figure 1. This is essentially a 
polarity counter with a null -zone and overflow. 

This machine, which can be called a "linear" machine, operates as 
follows. If any sample, sayX.^>d, then the state of the machine moves one 
step to the right. If X^< -d it moves one step to the left. Otherwise it re- 
mains in the same state. It remains in the end states unless, e. g. , X <-d 

K 

in state m. The machine initially starts in one of the middle two states with 
equal probability (assuming m is even and symmetrical statistics). If at 
time n, S^> — , Hj is decided. Otherwise H q is decided. This is a time- 
invariant machine since the transition rules from state to state do not 
change with k. Cover [ 1] has shown that, remarkably, as the length of the 
data sequence becomes infinite the probability of error Pe(m, n) asymptoti- 
cally approaches zero if m > 4 for time -varying machines. Indeed, this 
also holds true for m> 2 for certain statistics. Heilman and Cover [ 2] have 
found the asymptotic probability of error for time -invariant machines and 
have shown that this also goes to zero under certain conditions. They have 
also shown that the optimum machines resemble, asymptotically, that of 
Figure 1 but with d — oo . The convergence of probability of error is slow 
and little insight is gained from this asymptotic behavior about optimum or 
good machines and their performance for finite data sequences. 

Lynn and Boorstyn [ 3] have evaluated the performance of the "linear" 
machine for finite values of m and n. Typical results are summarized in 
Figures 2 and 3. Briefly these indicate that for ? bits of memory (m=8) 
performance close to optimum (infinite memory) can sometimes be achieved. 

These numerical results are for Gaussian statistics: f ^{x) Gaussian with 

2 ? 
mean p. >0 and variance or , f (x) Gaussian with mean -u and variance cr 

o 
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It is comforting that for this "linear” machine an explicit and simple 
expression for the probability of error can be obtained. This is given below 
in equation 2. 


( 2 ) 


F 1 " (m, n) = 


(a^-a j) 


1+ (a^ / a j) 


m/2 


m 


- 1 




j=i 

(j Odd) 


'-Pj 



aj =PrU k < -d| Hj } 

a 2 =Pr {|X k | <d| Hj } 
a 3 =Pr {X k > d| Hj } 

a l +a 2 + a 3 = 1 
and a 3 > a j * 

The first term in (2) is the steady state term since the second term van- 
ishes as n -* w . For a.j < a 2 << l (3jcan be tightly bounded by (3= 

1 XI 

for all j. Thus the second term can be bounded by -( a^-cij) (3 /(1-p) which 
is independent. of m. Thus m appears only in the first term and n in the 
second term. Both expressions are functions of a. and through them of d/cr 
and )jl / <r for Gaussian statistics. Thus (2) should yield considerable insight 
about the setting of the threshold d/cr and the dependence upon signal-to- 
noise ratio n/cr as well as the nature of the variation of Pe with m and n. 
Some results for a four state machine are given in Figure 4. 

To gain further insight we considered a machine with a different struc- 
ture which we call a majority rule machine. This looks only at the last 

2f-l non-null ( | | > d) samples and decides by majority rule. E, g. , if 

2£ - 1 

more than half of them are >d then is chosen. This machine has 2 
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states of memory and its performance is relatively easy to evaluate. 

This machine performs about as well as a linear machine with logarith- 
mically fewer (2 1) states but can serve as upper and lower bounds for the 
linear machine. (The lower bound may not always be applicable. ) Specifically 
if the superscript L denotes the linear machine and M, the majority rule 
machine, then 


A/f ? £ Hh 1 ^ L i* . » ^ . M i — 2 f ~ 1 v 

(3) Pe (2 ,n)<Pe (2*,n)<Pe (2 , n) . 

A comparison of these two machines is given in Figure 5. Neither of 
these machines are optimum, although Heilman haB shown that the linear 
machine for m=2 is optimum. 

Improved performance can be obtained if the machines are allowed to 
be time -varying . Muise and Boorstyn [4] have found the optimum time- 
vary-machine. It can be described as follows. If in state i at time k- 1 go 
to state j after receiving if 

MX.) 

(4) Vj (k)<«n f (x *j + 

L,(k- 1) is defined to be 

l 


(5) 


L.(k-l) 


. Pr{S k _,=i |Hj} 
= £n 11 

Pr{S k-l = l H o } 


and can be viewed as a likelihood ratio where the data is the only ob- 
servation that the machine can make, i. e. , its state, at time k- 1. Equa- 
tion (4) can be viewed as a time-varying quantizer where the Vj(k) are the 
threshold values and the L.(k-l) are the representation values. If, at 
time n, the machine is in a state S. for which L (n) > 0 then H is decided; 

otherwise H . 

o 

Algorithms have been obtained for calculating the sets of coefficients 
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L.(k-l) and ^(k) and for evaluating the performance of these optimum 
machines. The treshold values are given by 

f.(k) - f (k) 

(6) V j (k) = iTk)'- e (k) ’ 

j j 

where 

e j(k) = P {deciding Hj at time n | S^=j, Hj } 

f. (k) = P {deciding H. at time n I S. = j , H } 

J 1 1 k o 

It can be seen that if at time k the design of the machine in the past is 
known then the coefficients in (5) can be evaluated. However, the coeffi- 
cients in (6) depend upon the future design of the machine. The algorithm 
essentially starts with an initial design and iteratively runs back and forth 
in time applying equations (6) and (5) in turn. It has been proven that this 
algorithm converges and for typical examples the convergence is rapid. 
Figure 6 displays the coefficients for an example. Note that the coefficients 
for n=8 are not simply extensions of those for n=4 although the difference 
does not seem to be great. This suggests a suboptimum and simpler design 

which would find the best coefficients for each k as if it were the time of 
decision. 

Typical results for the optimum machine are also shown in Figures 2 
and 3 for Gaussian statistics. Also given there for reference is the per- 
formance of the optimum infinite memory machine and of the polarity coun- 
ter. These results show that remarkably good performance can be ob- 
tained with relatively simple machines and a very small memory- size 
limitation. 

The basic results in the optimization above can be inferred from the 
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following argument [ 5] . Let C denote a correct decision at time n. Then 
at time k 


1 m oo 

(7) P(C|S k _,=i) = L ^_ 1 / oo Pr(C l S k =£lS k-l =i ’ X k’ H j ) ' 

p r (S k ={ l S k-l =i,X k’ H j> • 
f Xk < x kl s k-i =i ' H j )Pr(H i |s k-r i)dx k 

But the first term in (7) is Pr(C |S^=5 , H^) and is used to form the terms 
in (6), the second term is Pr(S k =f |S k _ 1 =i. x^) = 0 or 1 depending on the de- 
cision rule at time k, the third term is I ^ 5 ) * an ^ ^ ourt ^ 

v J J 

term is used to form (5). If A u (k) is the region of ^ such that a transition 
will occur from S^_ ^ = i to S^= It, then 


m 


( 8 ) 


picis *-i 4 


(k) 


1 

V Pr(C |S. , H.)f.(x_ ) • 

u=6 k J J 

In 


Pr( H .|S k ,i = i) 


lx k 


The optimum transition rule for going from ^ to can be found 
from (8) by, for each i and x fe , finding the f that maximizes the bracketed 
expression. Muise and Boorstyn, in addition to other results, have given 
structure to this updating rule and proven its optimality. From (8) it is 
easy to extend this results to include time -varying statistics, many hypo- 
theses, and time-varying memory sizes. 

Although the memory has been constrained to be finite by limiting it 


to a sequence of states from a finite alphabet, the coefficients in (4) in the 
above analysis are not so constrained. (The arithmetic operation indicated 
in (4) is performed instantaneously, i. e. , not stored, so that no discrete 
limitation need be put upon it. ) These coefficients must be stored and there 
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are approximately 2nm of them! It is not known at this time what the effect 
would be if the storage of these coefficients were also limited. However the 
storage of the states is data sensitive (read-write memory) while the coef- 
ficients are (except for an adaptive machine) fixed (read-only memory). 
These differences can be exploited, especially by time sharing the coeffi- 
cients and the arithmetic operation with many similar detection processes. 

A sketch of such a structure is given in Figure 7. 

Further research is continuing on many aspects of this problem - op- 
timum time -invariant machines, algorithms for M-ary detection, other 
communications receivers, such as equalizers, more practical transition 
rule algorithms, time- sharing implementation, etc. It is hoped that these 
investigations will provide a new and more efficient viewpoint for designing 
communications systems incorporating fundamental digital constraints. 
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Figure 1. A " Linear" Machine 
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FIGURE 7: MACHINE CONFIGURATION - 
TIME SHARING' MODE 
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II. SIGNAL PARAMETER ESTIMATION FROM 
DISC RETE-TIME OBSERVATIONS 


This paper is intended to present some new results relevant to tlie prob- 
lem of estimating the parameters of several single -frequency tones from a 
finite number of noisy, discrete-time observations. This problem has ap- 
plication to data system testing, radar, and other measurement situations. 
The general model is indicated on Figure 1. In general the signal has 
k 

the form Y b, exp [j(w.t +0.)] . In a working system the imaginary part 
4J i i l 

may be derived from the real part by a Hilbert transformer or the imaginary 
part may not be processed at all. 

Time does not permit an examination of all the possible cases. There- 
fore, we will only discuss the case of a single tone whose imaginary part 
is a single cysoidal tone. An understanding of this case is fundamental to 

an understanding of the other cases. 

The parameter estimation problem was formulated and examined by 
Slepian in his often- quoted paper of 1954. ^ ^ Our intent here is to study, 
in more detail, a specific variation of the problem. 

The real part of the signal, s(t), is b Q cos(w 0 t+G 0 ). Suppose all three 
parameters are unknown, but only b^ and are to be estimated. (The 
phase can be estimated but we will not do so. ) 

The computer input will be two sets of samples, 


X = [ 


X 0 .Xi. 


X 


N-l 


] and Y = [Y 0 ,Y lf 


N 


-ii 


where 


X = s(t ) + W(t ), n = 0 to N - 1, 
n n n 

V v 

. y = B (t ) + W(t ), n = 0 to N - 1, 
n n n 


( 1 ) 


( 2 ) 
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and 

s(t) = b Q sin ( a) o t + e -0^ * ^ 

V 

W(t) is the Hilbert transform of the noise, W(t). We only consider the case 
of independent noise samples. 

If we write Z =_X + jY then the joint probability density function (p. d. f. ) 
of the elements of the sample vector Z when the parameter vector is a is 


given by 





ll = b cos (c*> t + 0 ) , (6) 

'll n 

and 

v = b sin (ut +0 ). (7) 

n n 

We will write t = (n n +n) T. 

n 0 

In developing the topic we will consider three main aspects of the 
problem. First we will examine lower bounds to estimation error, in 
particular Cramer -Rao lower bounds. Then we will develop and analyze 
maximum -likelihood (ML) estimators of the signal parameters. Finally, 
we will discuss practical estimation algorithms and simulation results. 
The frequency estimation algorithm has a threshold effect, which we will 
also discuss. 

BOUNDS 


Let us first look at lower bounds to estimation error. In an estimation 
(or measurement) system it is important to have numbers that indicate the 
best estimation that can be made with the available data (the observations). 



In a measurement system RMS errors are important and are often used 
as a measure of system inaccuracy. Estimation bias is of secondary im- 
portance, although it is generally desirable to minimize bias. In this report 
we will generally find that the bias can be neglected. Thus RMS errors will 
be the important consideration. 

Lower bounds to estimation accuracy have been studied by many people. 
Some of the better -known bounds are: the Cramer -Rao (C-R) bound t ^ and 
itB generalizations ^ 3 ^,the Bhattacharyya system of bounds ^,the Barankin 
bounds ^ ^,the Ziv-Zaki bounds ^,and the Chapman- Robbins bound 
The C-R bound is imbedded in the Bhattacharyya system of bounds and in the Bar auk 
in bounds . The Chapman- Robbins bound is related to the C-R bound. The difference 
is that the Chapman- Robbins bound avoids the need for the regularity conditions 
required in theC-R and Bhattacharyya bounds. Seidman has observed that many of 
the bounds are loose, especially at signal-to-noise ratios (SNR)above threshold ^ ^ . 
Some are also difficultto use. We will use maximum- likelihood{ML)estimation and 
will generally be able to keep the bias very small. Thus, above threshold, 
the C-R bound will apply. We will separately evaluate threshold effects. 

r 31 

The generalized C-R bound is due to Rao L J .It can be shown that this 
bound is the best (tightest) of a certain class of bounds ^ 

Before reviewing the bound it is helpful to make some definitions. We 
assume the indicated operations are legitimate. (They are in this problem. ) 

Let H(Z, a) = log f (Z_; a) - (9) 

Let £> be the vector defined by its typical element: 

a =5^-H(Z,a) . 

1 

The Fisher information matrix, J(a), is defined by 


( 10 ) 
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A typical element of J is 


J. . = E { H .H .} = - E { H : 
ij ai aj aiaj 

v v V y 


Let a ( Z ) be an estimator of a. . Let the matrix D be defined by typical 


element 


j 

Let the matrix C{_a) be defined by 


C(a) = E{a-a)(a-a) T } . (14) 

With these definitions and certain well-known regularity conditions on 

f(Zjki) and J, the generalized C-R bound is the statement that the matrix 
-IT 

C -DJ D is positive semidefinite. From this we find that 

Varf^} >[ DJ- 1 D T ]. i , i =1.2,3,. . . p. (15) 

where [ ] „ denotes the ith diagonal element of [ ] and p is the number 

of elements in a . We will ignore the conditions for equality to occur in 
(15) since they are generally not met in the problem at hand. 

E{ “j } = 6 ij ' (16) 

which happens when a is an unbiased estimator of a. , then D = I and the 
bounds become: ^ ^ 


Var { cL } , i = 1 to p . 


It can be shown that with f(Z;_a) as given above 


9a. 9a ; 
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where 


^0 = [“o' W • ' (19) 

After determining the elements of the J matrix, the details of which 
are omitted, one obtains: 


J = — 
• 2 

or 


b 2 T 2 (n 2 N + 2n 0 P + Q) 0 bJ'T^N+P) 


b y T(n Q N +F) 


N 


b i N 


( 20 ) 


where , 

N-l 

P- l 

n=0 


N- 1 


Q = n 2 = N(N^l)(2N-l) 


( 21 ) 


( 22 ) 


n=0 

and n Q T is the time at which the first sample is taken. 

From here on let us suppose n Q = 0. This will be convenient later. 

From (20) the bounds will depend upon n^ . We will not dwell upon 
this dependence. * 

22 

We can see from the zeros in J that the level bound, Var { a^} J * 

• , l i. _ a 1 . _ .1 ~ /(-»*• ^ t* Vtl HWn T f 

Will L>C tilt; b Ail it; wiicliicjl ui hvm. tnc puoot «■“«/«* — 

the phase were known the J matrix would be 


11 

0 


22 


for example. If the frequency were known the J matrix would be 



21 


J 


22 

0 


0 



Observe that the elements of J are not functions of u> or 0 . This is 
not true in the more general cases. 

Inversion of the J matrix is easy. When all of the parameters are un 
known the results are: 


Var {w A }> 


12tr‘ 


0 b 2 T 2 N(N 2 -1) 


Var { bg } > <r / N 


Var{ e 

U bJjN(N+l) 


( 23 ) 

( 24 ) 

( 25 ) 


MAXIMUM LIKELIHOOD ESTIMATION 


Now let us look at ML estimation of and b^, the frequency and 
level. The ML estimate of a. is the value of a, say a , that maximizes 
i{Z;a) when Z is the observed sample vector. 

The maximum of f(Z_;^) will occur at the maximum of log(f) =H( Z , a.) , 
or at the maximum of the function 


L = -H ' 


( 26 ) 


n 

Since T,X 2 , Y,Y 2 * and cr* 2 are not affected by o, b, or 0, we can drop them 
n n 


from L and use L^: 


r 

y,(X MY v ) - h J, H 2 + v 2 

n n l J 


( 27 ) 
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2b 

N 


Y f X cos(nwT+9) + Y sin(no>T+0)] - b' 
U L n n 

n 


Define 

F(Z.«)=Z (X n cos nwT +Y n sin nwT) 


n 


and 


G(Z,w)=^{Y cos nuT-X n sin nu)T) 
n 


Then 


Lj = [F(Z ,w) cos 0 + G (Z , co) sin 0] - b‘ 


(28) 


(29) 


(30) 


(31) 


Assuming b 0 > 0, we need to maximize F cos 6 + G sin 0 over both 
0 and <o . The maximum over 0 is >1 F +G . Thus the ML estimate of w q 
is the value of co that maximizes 


B(co) = [ F 2 (Z,<o) + G 2 (Z, «)] /N' 


The ML estimate of b^ is -*\ B( u) . 


(32) 


(31) 


Relationship to DFT 

Recall that the discrete Fourier transform (DFT) of the vector Z is 

, , . u [12,13,14] 

the set of complex numbers: L 

N- 1 

A 1 V z -j2*nk/N „ 1# - 2 N - 1 . (34) 

K N &0 n 


If we define a function A (u>) : 

N-l . m 

... 1 V ^ -jn<o T 

A(W) =N L. Z n e 
n=u 


(35) 


then 



2 3 


, k = 0, 1. . . , N- 1 . (36) 

It is easy to show that the function B(go) defined above is 

B(u) = |A(«)| 2 . (37) 

We see, therefore, that ML estimation of (Oq and is related to the DFT 
of the sample vector Z . This fact will suggest a practical algorithm later 

Properties of co ^ 

Because of the nonlinear nature of the estimation algorithm we are 
unable to derive the distribution functions for b Q and . We can, however 
show: 

1. The p. d.f. of go is symmetrical about modulo go 

0 s 

a 2 

2. Var {co } is proportional to co and independent of go (go =2*7 T. ) 

s Os' 

Noise Model . The following noise model leads to easy proofs. Let {V } be 

n 

a set of independent Rayleigh r. v. with parameter 1. That is 

I v exp [ -v 2 /2] , v > 0 J 

| (38) 

1 

0 , v< 0 ) 

Let be a set of independent r. v. uniformly distributed over 

(- 7 r) . 

We model the Gaussian samples as 



W 

n 


= ijV cos <b 
n ^n 


and 


W = a V sin <b 
n n T n 


(39) 


( 40 ) 
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Analysis . Recall 


Z = X + jY^ . 
n n n 


Then, using our noise model, 


jfnw.T + e,) j<l> 

Z = b.e 1 1 +cr V e n 

n 1 n 


(41) 


Thus 


V Jn " T 

n 


1“ 


N 


1 . — -j(nfl-<b +e, + nu ,T) i 

-[l b ie -^ +a V n e 1 1 ] 

n 


where (3 = (w-u) ^)T 


Let 


v = tb - G, - noo _ T , n = 0 to N - 1 . 
m T n 1 1 


(42) 

(43) 

(44) 


Since the <j> are independent and uniform on ( - 7r , tt ) , in effects so are the 


n 


■Yn • We are therefore justified in writing 


je, 


A(w) = s tt 2[ b o e ' jnp + ”‘ V n' 


-j( n P- yJ 




(45) 


Thus B(w)=-^ [ b 0 S cosn f3 +0 'E Vn cos ( n 0 - T n ) ] + •" 

N L n n 4 

, [ b 0 Z sin n P + o- Z V n sin ( n P“ V n )] 2 . (46) 


The independence from 0 q is now obvious. 

Without loss of generality, let (3 be the value of (3 in the range ( -tt, tt) 
that maximizes B(o>) . The estimator, w , will then have the value: 


w P 

w = oj „ + — _ s ■ — ; modulo w 

s 


0 0 27T 


(47) 
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Observe that B(w) is an even function of the pair p, y . <48) 

The statistics of - X are the same as the statistics of X ■ Thus the statists 
oi -x • Thus the statistics of -p must be the same as the statistics of p . 

Hence the p. d. f. of (3 must be an even function of p and E {p} - 0. From (46), 
the statistics of p do not depend upon or 6 q , but do depend upon the SNR, 

b 0 /2* 2 . 

Since we choose according to (47), the p. d. f. of p is related to the 
p. d. f. of u in the manner illustrated on Figure 2. The p.d. f. of u Q is even 
about u. except for the part from 2w to u g , when u Q < -^-(or part from 0 
to 2w - U , when u> > — ) . 


Consider the situation when w Q < -y- . If Pr{2w Q < w Q < w 8 } 18 Bma11, 
which it is when the SNR is large enough, the E{u } = u Q or u Q is unbiased. 
If Pr{2w < u < u s ) is significant then u Q is biased in the direction of 
„ s /2 . In other words, E{ u „-«* Q } > 0. If « „ > » ,/2 the above remarks 
reply with the obvious modifications. Observe that due to the symmetry of 
the problem, the bias of must be an odd function of u Q about u Jz. 

It is easy to show that if u Q is equal to zero, w Jz, »r« 8 the p. d. f. 

of ui - is even about w Jz. Thus, in these three cases E(io q ) = <■> g / 2 • We 
0 s 

see, therefore, that the bias of has the following values: 


W 


E(c> 0 -u> o> 


u> / 2 
s 


w / 2 
s 


w 


“ 2 /2 


Clearly we expect to make frequency estimation errors if o) Q is close 
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to zero or co . At moderate SNR, say above the threshold region, we would 

not expect trouble if the difference between and zero (or co^) is at least 

four times the RMS bound. In practice we do not expect the ambiguity at 

zero at co to be a problem, 
s 

The variance of (3 depends only upon the SNR. Thus the variance of co ^ 

2 . . 
is proportional to and is not a function of 0 q. The variance of co^ is a 

function of co^, but its variation with is small at SNR above the threshold 

region. Hence the C-R bound for unbiased estimators is appropriate in this 

region. 

An example. When there is no noise then it can be shown that 


A(«o) =b n e j9 °e- j(N - 1>Z 11 ^ ^- . 
0 N sin z 


(45) 


where 


<o -co 


z = 


-= T = 7T (co -w )/ co 

2 O s 


(46) 


B(co) is | A (co ) | . Thus 


? 

, v , 2 sin (Nz) 

B(u)) =b o~T S J ~ • 

N sm (z) 


(47) 


This function is shown on Figure 3. The function is symmetric about co 

7 

and has period co . The global maximum occurs at co and has value b 
s u u 

Notice numerous low-amplitude maxima. Without noise the ML estimates 
of co^ an b^ have no error. 

When noise is present B(co) loses its clean, symmetrical shape and the 
minor maxima get larger. The global maximum is usually close to co^. 
Figure 4 illustrates this situation. 

If the SNR is small B(co) will occasionally be so badly distorted that 



27 


the global maximum occurs at a frequency far removed from co Q as illus- 
trated in Figure 5. When this happens, the ML frequency estimation al- 
gorithm makes a large error. It is the occurence of these rare but large 
errors, which we call sports, at low SNR that causes VAR {« Q } to pull away 
from the C-R bound. 

AN ALGORITHM 

As indicated above, once an estimate of w Q is made, estimates of b Q 
and e Q can be done by straight-forward computations, using appropriate 
equations. Thus the difficult and time-consuming part of an algorithm is 
the part that locates the maximum of B(w) • This part is essentially a search 

routine. 

One way to develop an algorithm is to use a two-part search routine. 

The first part calculates B(u) for a set of w values between zero and to g , 
and identifies the w that maximizes B(«) over this set of to values. The 
second part locates the local maximum closest to the value of to picked out 
by the first part. We call the first part the coarse search and the second 
part the fine search . 


The Coarse Search 

For the coarse search we evaluate B(w) at the set of frequencies 
defined by 


_ 27Tk 
W k " MT 


k = 0, 1, 2 M - 1 , 


(48) 


where M is a power of 2 greater than or equal to the number of samples, 
N. We also always choose N to be a power of 2. Thus M/N is also a 

power of 2. 

Observe that the set {A(to k » is the DFT of the set [Zj defined by 
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Z n , n =0, 1,2,..., N- 1 

( 49 ) 

0, n=N, N + 1, . . . , M- 1 

The set {B(u> )} is simply the set { | A(co ) | } 
k h 

The output of the coarse search is the value of say that cor- 
responds to the largest member of the set {B{w )} . 

K 

It is natural to use M s N in the coarse search. However, it turned out 
that the w ^ thus obtained was the wrong choice (not close to the global max- 
imum) often enough at low SNR to cause trouble. We found that the number 
of wrong choices was significantly reduced when the coarse search used 
M/N equal to 2 or 4. 

Our coarse search uses a fast Fourier transform (FFT) algorithm to 
compute the desired DFT. 

The Fine Search 

The fine search algorithm locates the value of co closest to w that max- 

i 

imizes B(u). If the derivative of B(w) at B'(w^) is positive, the de- 
sired maximum is at a value of u greater than oo . Otherwise the desired 
maximum is at a value of co less than or equal to , 

Consider Figure 6. Given a frequency to start from the problem 
is to locate the closest zero in B'(oj) with B" (a) ) < 0. In the figure this oc- 
curs at point B. Points A and C correspond to minima in B(o>). 

Our fine search algorithm computes B'(u) at points P through P , 

18 

finally locating points on either side of the desired zero, points P and P . 

7 8 
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Then we use the s ecant method {Action [ 12] , p. 52) to compute successive 
approximations to point B and the frequency estimate, u . The iteration 
formula is simple: 


cj . 


i+1 


B'(cj. , )to . -B '(<*). )w . . 

l- 1 l i i-l 

-B'^) 


(50) 


The process is stopped if B ' (w = B 1 (cu , ^), if B'(w ) = 0, or if 

+i _0J i I * s ^- ess than 1 percent of the square root of the bound on the 
variance of w. . 

l 

If < u then the initial steps of the fine search work to the left of , 

using B' (co^) < 0 to indicate that w < . In either case the initial steps are 

to 

small enough to avoid missing the correct zero in B'(u >); the step size is 

5N 

Threshold Effect 

It is well known that nonlinear estimation is generally plagued by thresh- 
old effects. At low SNR there is usually a range of SNR in which the mean 
squared error (m. s. e. ) rises very rapidly as SNR devreases. The SNR at 
which this effect is first apparent is called the threshold. Receivers are 
often said to operate above or below threshold. 

Digital frequency estimation also has threshold effects, generally con- 
nected with the occurence of sports. In this section we>present a calcula- 
tion of threshold effects. The result accurately describes one particular 
model. 

Consider the estimation of the frequency of a single complex tone. Let 

the sampling frequency be u and the number of samples be N. Assume the 

s 

CO 

phase is unknown. Suppose the tone frequency is u = — — . Assume the 

« v Ct 

algorithm is the following, using M = N: 
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N-l 

1 v >7 -j2 TTnk /N . . 

1. compute the set: A^_ = /j Z^e . This i! 

n ~ 0 

the DFT of {Z } . 


n 


2. Identify the largest I A k I » sa Y l A jg I* and record £ . 

3. Start the fine search at w = m . and continue until 


the maximum of B(w) is found. 

Since we chose « 0 = X' A N/2 Sh ° Uld b * the largeSt - Tbat iB ’ tbe C ° ar9e 
search should give £ = N/2. If £ * N/2 we say a sport has occurred. If £ t 

N /2 the m. s. e. is greater than zero and less than the square of the distance 
toC °£+l * 

We will approximate the m. s. e, in this case by the C-R bound for an 

2 

unbiased estimator, which we designate * From equation (23), 


CJ 


2 

CR 


3w 


2tt 2 pN(N 2 -1 


(52) 


where 
P = 


0 / 2(3 


2 


(53) 


If a sport occurs, the outcome of the fine search will be any frequency 

between zero and o> . The p.d.f. is approximately uniform because the 

s 

signal has little influence. Thus we write the m. s. e. when a spo t occurs 
as 


2 

to 

sp 



(54) 


In the general case where the signal frequency is not equal to s/2, 


equation (54) would become 
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<0 


<o 


CO 


sp 


12 


to 


(55) 


The total m. s, e. is the weighted sum of the two contributions, 
m. s. e. = ( m. s. e. /sport) p + (m. s. e. /no sport)( 1-p) 


Let the total m. s. e. be to 


to 


“e “ P TT~ + (1_p) 


Then we have 

2 


3co 


s 


2 P ?r 2 * * * N(N 2 - 1 


(56) 


(57) 


The RMS error is 


to 


RMS 


= fi 


(58) 


Next we calculate the probability of a sport, p, and verify that when a 
sport occurs all possible t except the correct one are equally likely. 

Probability of a Sport . Let C = | A | , k = 0 to N - 1, (59) 

K iC 

where A^_ was defined above. When both signal and noise are present, 

to 

each is a random variable. If the signal frequency is ■ 8 and the 

2 

noise samples are independent, normal, and zero-mean with variance cr 

i it ca 

tioni 10 J 


then it can be shown that the C are independent with Rayleigh distribu- 

K 


W = la 


NCf' 

k 

NC. ’ 2 

k 2 cr 

e 


, c , > ° . 

k — 

k ± N/2 


( 60 ) 


Let r = N/2 . C has a Ricean distribution: 

r 


.[ 10 ] 


NC 

f (C ) =-r- 
r r 2 


N(C 2 +b 2 ) 

r 

•5 2 

Zcr 


NbC 


, C > 0 
r — 


( 61 ) 
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where b is the signal amplitude. Then 

1 - p = P { all C < C } 
r r k r 


= f F {all C < C /C = x } P {C = x} dx . 
Jr k r r r r 

x 

But P {all C < C fC = x} = [ P {C < C /C = x l] 1 *' 1 
r k rr L r 1 rr J 

Thus 


1 - P = 


fo 


f r (x) 



x x __ 2 2 , 

/ f <y)dy = / e /2r 

0 K 0 cr 


= 1 - e 


-Nx 2 /2t r 2 . 


(62) 


(63) 


(64) 


Then 


1 - p = 



r 2 “i 

NX 

N- 1 

N{X 2 +b 2 



„ 2 


2 

► « 

f 

. 2c r 

1-e 

NX 
_ e 

2 cr j 

NbX 

J o 

L 

2 

<r 

o 

2 

cr 


dx. (65) 


After some further work, we obtain 


N-l 

i-p = Z 


k=0 


(N-l) 1 (-l) k -Npk/k+1 
(N-l-k) ! k! (k+1) 


and 


N 

Z 


N!(-l) 


m 


- Np 


m- 1 
m 


m 


_ 2 {(N-m) 1 m! 


( 66 ) 


(67) 


It is easy to verify that the limit of p as p — *0 is 1 , which is in agree- 

ment with our assumptions. The formulas for p given above are neat 
closed forms that cannot be summed on a computer because the terms 
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Nt(-l) 


m 


(N m)!m! ver ^ l ar 8 e and- alternate in sign. Thus to compute p it was 

necessary to use the integral form and do numerical integration. When Np 


is large one can use the first term of p: 


„ _ N-l -Np/2 „ < in -3 

p e , p < 10 


( 68 ) 


The calculated values of p are shown on Figure 7 . 

Approximate RMS Frequency Error . We uaed the above formula for several 
values of N as shown on Figure 8. The small circles of the curves represent 
the results of simulations. As can be seen, the simulation results agree 
with the calculated curves. The curves are similar to the well known results 
for the continuous observation case. See Van Trees, p. 285. 

On Figure 8 the simulation for N = 16 each were done with 800 estima- 
tions. The ones for N = 128 were done with 500 estimation!. Five hundred 
estimations will have a sample variance within 12. 5 percent of the true var- 
iance with the same confidence. 

One would not usually operate a system at SNR below the threshold. 

Thus Figure 8 is useful mainly because it shows the SNR at which the thresh- 
old effect starts. All SNR above threshold can be considered to be 'high 
SNR U in the sen.se that the variance of ML estimators equals the C-R bounds 
at high SNR . 

Level Estimates 

The simulations described above included level estimates according 
to equation (33). In every case the RMS level errors were almost equal 
to the C-R bounds. Threshold effects were not observed! 

Remarks 

We ran the above -described simulations with M = 4N instead of M= N. 
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There was no significant difference in the results. Since using M= 4N is 
more likely to result in correctly locating the global maximum in B(co) , we 
are led to believe that Figure 8 truly depicts ML estimation when the fre- 
quency is one half the sampling frequency 

The next question is, what about different signal frequencies? We ran 
the simulation with M = 64, N = 16, and f Q = 2120 Hz, using -10, -5, 0, and 
5 db SNR. The only point different from the 0 points is the Q point on Fig- 
ure 8. As before, level estimates did not show a threshold effect. 

SUMMARY 

This has been a quick trip through a study of the problem of estimating 
the frequency and level of a cysoidal signal from a finite number of noisy 
observations of the signal. We derived the equations that describe the Cramer - 
Rao lower bounds to the variance of estimation errors. Then we derived 
the maximum -likelihood estimators and showed their relationship to the dis- 
crete Fourier transform. The analysis of the ML estimators revealed some 
of their properties. Then we looked at an algorithm suitable for implementa- 
tion on a digital computer. The algorithm almost always yields ML esti- 
mates. We were able to derive an expression for the threshold behavior 
of the algorithm. Simulation results verified the analysis. 

The overall conclusion for the case studied is that ML estimates is 
feasible and will yield estimates which are as good as permitted by the C-R 
bounds (above threshold). 

The general cases of real tones (sinusoidal signals) and of many tones 

are, in a sense, an extension of the case studied here. The presence of 

several cysoidal signals introduces complexity in the bounds, ML estima- 

[ 10 ] 

tion, and practical algorithms. These matters have been studied 
are not reported here. 


but 
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FIGURE 1 - PARAMETER ESTIMATION SYSTEM MODEL 
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FIGURE 2 - Relationship of p. d. f . of j-5 to p. d, f. of u . 
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FIGURE 3 - Shape of B(w) from single complex tone without noise. N is 16 
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PR03, OF A SPORT 

F 1 ,* 1 • I l . 
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FIGURE 8 - Approximate perfor*»nce of ML frequency estimate of single 
complex tone at 2004 Me* with unknown phase. 
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HI. DIGITAL FILTERING FOR RADAR SIGNAL 
PROCESSING APPLICATIONS 

A novel approach in synthesizing digital filters for signal processing 
applications is presented. This approach is an extension of the frequency 
sampling method of nonrecursive filter synthesis. Appropriate time delays 
are used in conjunction with a set of parallel complex exponential resonators 
whose outputs are summed to yield a desired filter impulse response. This 
synthesis method takes advantage of the known signal waveform structure 
and results in many fewer digital computations as compared to convolution 
processing. This approach is particularly suited to synthesis of matched 
filters for radar signal processing and yields matched or approximately 
matched filters which simultaneously have very low storage and very low 
computational requirements, 

I, Introduction 

Interest in the application of digital filtering to signal processing is 
increasing and numerous publications (Refs. 1-5) have been written on various 
approaches to digital signal processing. In this paper a novel approach to 
signal processing digital filter synthesis is presented which takes advantage 
of the signal waveform structure and results in a reduction in the required 
digital computations. 

The signal to be extracted from an interference background often con- 
sists of only a few data points as compared to the amount of data to be pro- 
cessed. An example of this is found in radar signal processing where the 
data to be processed usually consists of thousands of samples, whereas 
the signal is composed of several hundred or fewer samples. 

The input data to a radar analog signal processor consists of analog 
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voltage variations representing signal plus noise {or noise alone) corre- 
sponding to a given time or range extent. To use this information in a digital 
processor, these voltage variations should be sampled such that there is at 
least one sample per range resolution cell (defined as the reciprocal of the 
signal bandwidth). These samples are then converted into binary words which 
represent the amplitude and phase, or in-phase and quadrature components 
of the sampled analog voltages. These complex valued words or samples 
are then operated upon to extract the wanted signal from the background 
noise. 

The radar signal to be extracted from noise generally has a time extent 
less than that of the input data of interest, but longer than a range resolu- 
tion cell if the signal bandwidth -time (BT) product is greater than one. The 
presence and location of the signal in the input data can be determined in an 
optimum manner by convolving the complex conjugate of the sampled signal 

reversed in time, s*(-n), with the input data,y(n), as follows"^: 
n 

a(n) = T, x<J) s * [N-n+j] ; n = 0, 1, 2, . . . , J- 1 (I) 

j=0 

where a(n) is the sampled cross-correlation function of y(n) and s*(-n). It 
is assumed that there are J range data samples and N signal samples. If 
the signal is present at some point in the input data, the cross -correlation 
function^ will peak up at that range location. 

^It can be shown that, for purposes of signal detection and signal parameter 
estimation, the signal can he extracted from Gaus sian-distributed inter- 
ference in an optimum manner by matched or correlation filtering (Ref. 6). 
This type of filtering involves the convolution of the signal and data as 
described above. 

^Since the y(j); j 1, . . . ,n in (l)are random variables, then each afnl is 
also a random variable. Therefore, the discrete function {a(n);n=0, 1, 

...» J- 1} is a sample from an ensemble of discrete random functions. 
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This type of signal processing is referred to as ’’Moving Window Con- 
volution Processing” since the finite -length samples signal reversed in time, 
s*(-n), acts as a window moving past the range data. Note that each value of 
a(n) is formed by one position of the ’’moving window" relative to the range 
data. If J samples of range data are present, there must be J different 
positions of the”moving window" in order to produce J different values of 
a(n). 

Another method of extracting the signal from the range data is to per- 
form the filtering in the frequency domain. That is, transform all. ot the 
range (time) data to the frequency domain by use of the Digital Fourier 

Transform (DFT) as follows: 

J-l 

X(k) = Y. y(n) exp f -j 1 , k = 0, 1, 2, . . . , J- 1 (2) 

n=0 l j j 

Then, perform the filtering by multiplying the freqpency data by the 
Fourier Transform of s*(-n): 

Y (k) = S^(k) X(k); k = 0, 1, 2, . . . , J- 1 0) 


The resulting frequency samples are then transformed back into the time 
domain by use of the Inverse Digital Fourier Transform (IDFT): 


a(n) - 


, J-l 

1 Y 
r Li 
k=0 


Y(l<) exp j 


7rkn 




n = 0 , 1 , 2 , 


.T - 1 


(4) 


It the signal is present, the c ross -cor relation function 1 a(n), will peak 
up at that range location. 

If this method of processing were used, the Fast Fourier Transform 
(FFT) and Inverse Fast Fourier Transform (IF FT) would probably be used 
instead of the DFT and IDFT as indicated above. The transformations 
accomplished by the DFT and FFT, or the [DFT and IFFT, are identical; 


|| See note on previous page. 
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the way in which the transformations are accomplished is different (Ref. 7). 

This type of processing is referred to as "batch processing", since the 
complete set or "batch" of range data is processed at once rather than in 
subsets as in moving window convolution processing. 

Batch processing is sometimes more desirable than moving window con- 
volution processing because many fewer digital computations are required 
when the Fast Fourier Transform (FFT) is used. However, the disadvantage 
of batch processing is that a large amount of storage is required to store the 
thousands of range data points. To date, batch processing for many radar 
applications is too expensive to implement because of this large storage 
requirement. 

On the other hand, moving window convolution processing requires less 
storage (only the number of signal samples) but more digital computations 
and is, therefore, also very expensive to implement. However, if the 
number of required computations in moving window processing could be re- 
duced significantly, then the combination of relatively low storage and low 
computational requirements would make it a desirable approach to digital 
signal processing. Hence, it is of interest to consider ways of reducing the 
number of digital computations required in moving window processing. 

In the following sections of this paper, a synthesis technique is dis- 
cussed which takes advantage of the signal waveform structure such that 
the number of digital computations required is at least an order of magni- 
tude lower than that required in moving window convolution processing. 

In addition, the storage required with this technique is about the same as 
in moving window convolution processing. 
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In order to explain this technique clearly, several methods of synthe- 
sizing digital filters are briefly , reviewed in the next section of this paper. 

II. Methods of Nonrecursive Filter Synthesis 

In the digital filter synthesis problem, the desired impulse response 
is assumed to be known and is represented by the sequence 

{ h(n): n = 0, 1, 2, . . - } (5) 

where each h(-) is complex valued. This discussion will be restricted to 
the synthesis of nonrecursive digital filters which are characterized by a 
finite -duration impulse response such as is required in-moving window pro- 
cessing. Hence, the impulse response will be represented by a sequence of 
N complex numbers 

{ h(n): n = 0, 1, 2, ...» N-l} (6) 

A straightforward method of synthesizing a filter with this type of im- 
pulse response is the convolution or "tapped delay-line" filter realization 
illustrated in Figure 1. In this realization the N weights are the complex 
values of the N samples in the filter’s impulse response, and each of the 
N-l delays used corresponds to the delay between the filter impulse re- 
sponse samples. 

Now, if the impulse response is required to be 
h(n) = s*(N-n); n s 0, 1,2, , N-l (7) 

then the output of the convolution filter of Figure 1 is related to the input 
data, y(n), as gi ven hy Equation (1). Hence, this filter would perform the 
moving window processing desired. 

* In radar /sonar signal processing the number sequences encountered 
are complex valued, hence, the need to assume h(-) is complex. 
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The number of complex multiplications required with this filter is N, K 
where N is the number of samples in the signal (and in the filter impulse 
response), K is the number of data samples to be processed, and where the 
number of data samples is assumed to be greater than the number of signal 
samples (K>N). 

Now let us consider another method of synthesizing a digital filter with 

the impulse response of Equation (6). 

Since the impulse response is of finite duration, it can be represented 
in terms of its Discrete Fourier Transform (DFT) { H k :k =0, 1, 2, . . . ,N-l} 
as follows: 


k(n) = ^ Y H k exp [ ^ ^ N^ n 1 ' n * 0, 1, 2, , N- 1 

n=0 


( 8 ) 


The z-transform of the sequence of (6) is 
N-l 

H(z) = Yj 1*(n) z n 
n=0 


(9) 


Substituting (8) into (9) and interchanging sums, the sum over the n index 
can be evaluated in closed form so that 

-N 


»•> H 

h, *”L — 


1 - z 


(10) 


1-z" 1 e*p[ J ^|f L ] 

Evaluating (10) on the unit circle where z = exp[j27rf-r] leads to the 
frequency response 

N-l T t 

sin Ntff 


N-l 

H(f) = exp [ - j(N- l)7rf t], Y n 

k=0 


' ^k . r -.i^rk 1 

ex Pl N J* 


iin 5r[ ff- | ] 


(U) 


The symbol t denotes the time spacing of the number sequence in (6) and 
is often assumed to be unity. Note that when z = exp | ^ ”■ ]> H(z) - 
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or, equivalently, when f = ^ ^ , H(f) = . 

The interpretation of (10) and (11) is as follows: the frequency response 
of the digital filter which is described by the complex impulse response in 
(6) is a sum of frequency responses, each with a complex weight 

k 

H k and a center frequency — t where k = 0, 1, 2, . . . , N- 1, as illustrated in 
Figure 2. A realization of the filter of Equation (10) is shown in Figure 3. 
Each block in this figure represents a zonal filter with the appropriate cen- 
ter frequency and bandwidth. Since the term (l-z“ N ) is common to all of 
the filter functions, it can precede the filter bank and be shared. 

The parameters needed to completely specify the digital filter are either 
the sequence { k = 0, 1, 2, . . . , N~ 1} and the center frequencies 
^ £ k = "NT ' 1 = °» 2 » * • * » N- 1} or the sequence {h(n);n = 0, 1, 2, . . . , N- l] 

and the time delay between samples, T , If the desired frequency response 
is known, the filter weights are determined by sampling the desired frequency 
response at the frequencies f ; the complex values of these frequency samples 
are then used as the filter weights, On the other hand, if the desired 

impulse response is known, the filter weights, H f can be determined by use 
of the following equation, which is obtained from (8): 

N-l 

H k = Z H < n > exp [^fp 2 ] ; k = 0, 1, 2 N-l (12) 

n=0 

Note that the impulse response of the filter of Figure 3 is exactly that 
given in Equation (6) if the. complex weights are selected correctly. In other 
words, the samples taken from the frequency spectrum of h(n) are adequate 
to completely describe the filter frequency spectrum. 

It will be useful for later discussions to consider in more detail the 
impulse response of the realization given in Figure 3 . The impulse re- 
sponse of the £ block or subfilter of this filter is 
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exp[ J 2 ^ -- ? -]; n= 0, 1,2, . ■ ■ ; < = 0, 1, 2 N- 1 (13) 

-N 

which fs of infinite -time duration. However, when the term (1-z ) is in- 

til 

eluded with the £ subfilter such that the z -transform of that filter is 
, -N 

— (H) 

l - z lex p[^KT“] 

the corresponding impulse response is 

exp[j2^£.];n = 0, 1, 2, .... N-l ;< = 0, 1, 2 N- 1 (15) 

which is of finite-time duration. Therefore, the impulse response of each 
subfilter or block of Figure 3 can be thought of as a time response with a 
rectangular amplitude and linear phase function that is sampled at the time 
instants n T , n = 0, 1,2,..., N-l. The overall impulse response of the 
digital filter is the weighted sum of the N time sequences described by (15), 
where £ runs from 0 to N- 1. The composition of the composite filter re- 
sponse is illustrated in Figure 4, where N is taken as four. 

In some applications, one would expect some of the complex weights 
{ H ; k =0, 1, 2, . . . , N - 1} to be zero. For example, if the desired impulse 
response {h(n);n =0, 1, 2, ... , N- 1} represented data with a small band- 
width (the complex value of the sequence varies slowly with increasing n), 
many of the complex weights would be zero. A very special case would be 
where the magnitude of the complex impulse response did not vary with in- 
creasing n and the phase of the impulse response increased with n in a 
linear fashion; that is, the desired impulse response is given by 

H k ‘ expp£^];n=0, 1. 2, .... N- 1 (16) 

where the index k may take on any one of the discrete values 0, 1,2, . . . , N-l. 
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The resulting filter realization would have only one block or subfilter, as 
indicated in Figure 5. 

The number of complex multiplications required in the filter realization 
of Figure 3 is 2 • N • K, where N is the number of signal samples (and filter 
impulse response samples), 1C is the number of data samples to be processed, 
and K is asumed to be larger than N. The number of complex multiplications 
required for each subfilter (including the complex weights) is 2* K. There- 
fore, assuming that all N of the eubfilters are needed, N? (2. K) complex 
multiplications are required in the composite filter realization. 

In comparing the two filter realizations discussed thus far, one concludes 
that the use of the convolution filter realization of Figure 1 requires fewer 
multiplications than the frequency response realization and is, therefore, 
simpler and cheaper to implement. However, a further consideration of the 
frequency response filter realization leads one to a third filter realization 
that requires many fewer multiplications than even the convolution filter for 
some types of radar signals. This is the subject of the remainder of this 
paper. 


III. Synthesis Of Specialized Digital Filter Response 

In light of the discussion of the previous section, consider the digital 
filter illustrated in Figure 6. This filter is very much like the filter of 
Figure 3 except for the time delay following each subfilter. However, the 
impulse response of this filter is much different from the filter of the pre- 
vious discussion. 

The z-transform representing this filter function is 


F(z) 


Li-l 


-X 

I =0 


-£N/L . F f 

! ¥ 
L 



J 


(17) 
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where N is the number of samples in the filter's impulse response, L is the 
number of subfilter, and is assumed to be an integer. Evaluating ( 17) on 
the unit circle where a = exp (j2*fT) leads to the frequency response 


F(f) = exp -I* 


m 


-i 


L.-1 


i T 


, y £ 

[l ] 


exp | -jTT 


2£f T 


[l) 


it 

N 


sin [z] 


TTfT 


U8) 


sin 7T £ f T - 


t h 1 

N J 


When z = exp 


j2*< [-&] 


, F(z)=F.;or, equivalently, when f = [ , F(f) = F { 


The inverse z-transform of (17) yields the impulse response of the 
digital filter in Figure 6. 

1< F 0 ' F ,!n= 0,1,2 x- 1 ! 

i27rnL( % 

2N 


< F r w )-e 


]2nnL 

n . n n JL'+i .. 

* n L ’ L ' ’ 


- 1 


i4?rnL 


. L, N 2N 2N t 1 . l 

TT) e ;n = -T— , -TT +l# L 


(F 


2 N 


k*- «-•** H 11,1 


The magnitude and phase of this number sequence is illustrated in 
Figure 7. Those familiar with waveforms used in modern radar will re- 
cognize this filter response as that which is required to optimally process 
a Stepped Frequency Modulated (SFM) waveform (Ref. 8). As will be illustrated 
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in the following example, the filter of Figure 6 optimally processes, in a 
moving window fashion, the SFM waveform with many fewer multiplications 
than required by the convolution filter of Figure 1. Hence, Equation (17), 
which is realized as the filter of Figure 6, represents a method of filter 
synthesis which, for the SFM waveform, fulfills the goal set forth in the 
Introduction : to reduce the number of multiplications required in moving 
window processing. To illustrate this concept, consider the following ex- 
ample. 

IV. Optimal Filtering of the Stepped Frequency Modulated Waveform 

Suppose it is desired to synthsize a digital filter that is matched (in the 
signal processing sense) to the stepped frequency modulated waveform of 
Figure 8. In this figure, N and L are assumed to be 1024 and 16, respectively. 
The bandwidth of this waveform is ■ , where T is the basic time delay of 

the filter. The magnitude of the frequency spectrum of this waveform is 
shown in Figure 9, and is repetitive with period^ . The phase variation 
with frequency, which is not shown, is nearly quadratic. 

The filter which is matched to this waveform is shown in Figure 10. 

Its impulse response is the conjugated time-inverse of the waveform shown 
in Figure 8. The filter f s frequency response is the conjugate of the wave- 
form spectrum shown in Figure 9, and is composed of 16 sin Nx/sin x 
responses, C3-cli with nnitv Fi cin rp 11 i 1 1 n p c +■ 1 ^^ ^ ai*s 

of this filter spectrum. The ripples in the pass band and stopband are due 
to the sidelobes of the various sin Nx/sin x responses. Also note that nulls 
exist in the stopband of the filter where every sin Nx/sin x response has 
zero value. 

The time response of this filter to the stepped frequency modulated 
(SFM) waveform of Figure 8 is shown in Figure 12. This response, which 
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is the convolution of the impulse response of the filter of Figure 10 and the 
SFM waveform illustrated in Figure 8, is very similar to a sin x/x. function 
in the vicinity of the mainlobe response with the accompanying high sidelobes 
If it is desired to suppress these sidelobes, than amplitude weighting can be 
used. That is, the filter weights F ^ ; f = 0, 1, 2, ...» 1 can be selected so 

that the time response is similar to that in Figure 12, but with lower side- 
lobes. As an example, suppose the sin x/x sidelobes are to be lowered to 
30 dB below the mainlobe. Sixteen weights which are obtained by sampling 
the Taylor Weight function (Ref. 9) can be used in the filter of Figure 10 to 
yield the time response of Figure 13. 

The number of multiplications needed to process a sequence of samples 
with the filter of Figure 10 is L- K, or 16K where K is the number of samples 
to be processed. If the filter realization of Figure 3 had been used, the nufn 
ber of multiplications required would be 2 * N • K, or 2048K-. Hence, the 
savings in multiplication for this particular example are a factor of 128. If. 
nonunity weights are used with the filter of Figure 1Q> the number of multi- 
plications required will be 2 ■ L : K and the savings factor will be 64. 

Of course, the filter matched to the waveform of Figure 8 could have 
been represented with fewer than 1024 points, since the waveform is over- 
sampled by a factor of four. For example, the filter impulse response 
could have been represented by 256 points, in which case the resulting 
multiplication savings factor would be 32. However, the frequency response 
would look different than that shown in Figure 9. The waveform spectrum 

J 

would be spread out over the repetition period so that the filter would be 
essentially an all-pass filter with an approximately quadratic phase function, 
and would have to be preceded with a linear-phase low-pass filter. 
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V. Optimal Filtering of the Linearly Frequency Modulated Waveform 

This technique can also be used to synthesize digital filters matched to 
other important waveforms. For example, of great importance to radar 
signal processing is the linearly frequency modulated (LFM) waveform 
(Ref. 8). The magnitude and phase characteristics of this waveform are very 
similar to those of the stepped frequency waveform illustrated in Figure 8. 
The big difference in the two is that the LFM waveform has a continuously 
quadratic phase function in time, whereas the SFM waveform has a piece- 
wise continuous phase function that is only approximately quadratic. Be- 
cause of the similarities of these two waveforms, one would expect that the 
filter of Figure 6 could be modified to optimally filter the LFM waveform. 

If in the filter of Figure 6 there were more than one subfilter preceding 
each delay, as shown in Figure 14, than greater freedom in selecting the 
complex-valued impulse response would exist. In fact, if there were N /L 
subfilters preceding each delay, then any complex impulse response defined 
by (6) could be achieved by selecting the appropriate weights. 


| F { k ;f= 0, 1,2 L-l; k = 0, 1, 2 [£-] - 1 j 


( 20 ) 


which are related to the impulse response 
follows : 


h(n) ; n =0, 1, 2, . . . , N-l } as 
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i ,k 


U+l) 

- I 


n=* • 
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- 1 


h(n) exp 


- i 2 Trim 
N/L 




< = 0,1,2 ,L-1 (21) 


k =0, 1, 2, . . 



-1 


Hence, specific weights can be chosen so that the filter of Figure 14 is 
matched to theLFM waveform. 

The number of multiplications required to process a sequence of input 
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samples with the filter of Figure 14 is K- |^-j— J per subfilter for a total of 
2. K- N for the composite filter (including the N frequency weights). This 
is the same number of multiplications as required by the convolution reali- 
zation of Figure 3. Consequently, this realization has no advantage over 
that of Figure 3 if a filter perfectly matched to an LFM waveform is required. 
However, if a filter that is approximately matched to the LFM waveform is 
acceptable, then this synthesis approach is advantageous. For example, the 
filter of Figure 10 which only requires 16* K multiplications is almost matched 
to an LFM waveform of length 1024 T and bandwidth l/(4T). The phase func- 
tion of the SFM filter is not exactly quadratic, as is desired for a filter 
matched to an LFM waveform. However, the mismatch is small, and for 
some applications the SFM filter of Figure 10 might ge adequate for an LFM 
waveform. 

If a filter which is sligh'ly mismatched to an LFM waveform is accept- 
able, but less mismatch is desired than that provided by the SFM filter, 
then the filter of Figure 15 might be considered. In this filter, three sub- 
filters per subsection are used. The complex weights can be chosen to ap- 
proximate the LFM quadratic phase function, but with less mismatch than 
with one subfilter per subsection as in Figure 10. The number of multipli- 
cations needed for this filter is 3 • L • K, or 48 ♦ K. Hence, the savings in 
multiplications for this filter, as compared to the filter of Figure 7, are 
2048/48, or approximately 43. 

It is obvious that the LFM quadratic phase function can be approximated 
with increasing accuracy by increasing the number of subfilters per subsec- 
tion. Of course, the multiplications savings factor decreases with the in- 
creasing number of subfilters. The limit of this approach is when all (N /L) 
subfilters per subsection are used, in which case the filter is matched to 
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the LFM waveform and 2-N*K multiplications are required. 

VI. Conclusions 

Nonrecursive digital filters can be synthesized by several techniques, 
as described by various authors. Furthermore, the complexity of the filter 
may depend on the synthesis technique chosen. This paper has presented a 
technique for designing some types of finite -duration impulse-response di- 
gital filters which require many fewer multiplications than other known 
techniques. Any desired finite -duration impulse response can be synthe- 
sized with this approach. However, the savings in multiplications depend 
on the impulse responses desired. For some well-known radar applications, 
such as pulse compression, the multiplication savings factor can be very 
large. The examples discussed in this paper indicate that the savings in the 
required multiplication rate can be factors of a hundred or more. 
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Figure 1. Convolution Realization of Nonrecursive Digital Filter 
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Figure 4. Composition of Digital Filter Impulse Response (N = 4) 
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Figure 6. A Specialized Digital Filter Realization 
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Figure 8. A Sampled Stepped Frequency Modulated Waveform (N - 1024) 
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Figure 11. Composition of the Frequency Spectrum of Figure 9 
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Figure 13. Time Response of Digital Filter to Step Frequency Modulated Waveform 
Taylor Amplitude Weighting (-30 dB f n = 6) 




Figure 14. A Digital Filter Realization with L Subsections 
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Figure 15. Digital Filter Approximately Matched to the Linearly Frequency 
Modulated Waveform 
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IV, ON THE OUTPUT OF MULTIPLE SERVER QUEUES 
WHERE NOT ALL SERVERS ARE IDENTICAL 

It is known that if the input to an N server queue, with the servers 
having exponential service time distributions with identical means l/o-, is 
Poisson with paramter X, then, after attaining equilibrium, the output is 
also poisson with parameter X . It is shown in this paper that as long as 
the server a have exponential service time distributions, even if not iden- 
tical means, and equilibrium is attained with a Poisson input of parameter 
X , the output will also be Poisson with parameter X . Furthermore, for 
the output of each individual server to be Poisson, it is necessary and 
sufficient that the propability of a particular server being busy, conditioned 
on the number of customers in the system, be equal to the probability of 
any other server being busy, under the same conditioning. The output from 
server i will then have parameter X ■— - and be independent of the output 
from the other servers. However, for such equiprobable conditions to 
exist, restrictions must be placed on the dispersion of the mean service 
times of servers. 

I. Problem Description 

The problem of describing the output of a queue has been given con- 
siderable attention, primarily due to the problem's application to tan- 
dem queues. However, in most of the literature that attention has been 
focused on queues where the N servers have identical service time dis- 
tributions and on the output of the total system rather than on the output 
of any individual server. In this paper the output of an individual server 
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is examined. The motivation for this examination is the consideration of 
tandem queues in networks where the output of each server may feed into a 
different queue than the queue fed by the output of another server. In parti- 
cular, modeling of communication networks for message transmission from 
a terminal to a set of processors over channels of fixed capacity may be 
posed as an arrangement of tandem queues where the channel capacities and 
message lengths determine the various service times of the servers. 

The general queueing system for this paper is shown in Figure 1. The 
input is Poisson with parameter X. Each server i has an exponential service 
time distribution with parameter o\ , not necessarily all the same. P. J. 

Burke ^ has shown that if the service time distributions are all identical and 
N 

Y, cr. = No- > X 
i=i 1 

and equilibrium has been attained, the output of the overall system will be 

2 

Poisson with parameter X, same as the input. Edgar Reich has also shown 
this result by arguments different than those of Burke. In this paper the 
basic strategy of Burke is followed to show that the output from the overall 
system is Poisson even if the servers are not identically distributed, so 
long as equilibrium is attained and the system is Markovian. Furthermore, 
for the output of each individual server to be Poisson, it is necessary and 
sufficient that the probability of a server being busy, conditioned on the 
number of customers in the system, be equal to the probability of any other 
server being busy, under the same conditioning. The output from a parti- 
cular server i will then have parameter X ~i and be independent of the out- 
put from the other servers. However, for such equiprobable conditions to 
exist, restrictions must be placed on the dispersion of the mean service 


times. 
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XI. Equilibrium Probabilities 

In order to pursue the analysis it is first .necessary to determine the 
equilibrium probabilities describing the queue. The probability of a parti- 
cular server j being busy conditioned on i customers in the system will be 
represented as 


Q * = Pr (customer in cr . | 
j y 


i customers in system) 


The m,ean service time of the system, conditioned on i customers in the sys 

tern, is . defined as 
N 

j=l 


Finally, define the equilibrium probability of the state of the system as 
= Pr ( i customers in system) 

The equilibrium relations describing a system of N servers may now be 
written as below. 

XP 0 =>, Pj 


<X + H i )P l =*P ul + K i+1 P i+1 l<i<N-2 


(1) ( ;v+ M- n .i) p n .i = * P N-2 + * P N 


(X+o-) P. = XF. . + <r P. , 
J J' 1 J+l 


N < j < oo 


These equations are subject to the constraint 


00 


(2) Z P. = 1 
i=0 1 
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Define C = — , £ . - — , 1 < i < N- 1 , and = 1 . The solution to the above 
o- l . — — 0 

equations (1) subject to the constraint (2) is then easily found to be 

j 

p. = n g. p n l < j < n- l 

j . A i o — — 

J 1=0 

. 4.1 M N-l 

(3) P. = (e) J n 5.P 0 N< j< CO 

■* i=Q 


p o = 


ijJL 


N-l 


N-2 , 


2 ? «i] -q l, Kh 


j=o uo 


1=0 i=0 


III. Output of General Queue 

The equilibrium probabilities may be interpreted as the probabilities 
of finding the system in a particular state at any randomly selected instant 
of time. Central to the developments in this paper is the result that the 
probabilities of finding the system in a particular state at a randomly se- 
lected instant of time drawn from the set of instants immediately after de- 
partures are also equilibrium probabilities. This result is stated below as 
a lemma. 

Lemma 1 

The state of the general queue, subject to Poisson input and having 
attained equilibrium, after a departure is described by the equilibrium 
probabilities. 

With replacing T for 1 <_ i 5, N - 1 the proof of the lemma is completely 
analogous to the demonstration of a similar lemma by Burke. 

The above lemma and replacements can be used to demonstrate the 
following theorem which is a simple extension of Burke's theorem. 
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Theorem 1 

The output of the general queue with Poisson input having parameter X, 
and after attaining equilibrium, is also Poisson with parameter X . 

The demonstration is again completely analoguous to the demonstration 
by Burke. 


IV. Output of Individual Server 

The principal theorem of this paper pertains to the output of the indi- 
vidual servers in an N server queue. The theorem is stated below. 


Theorem 2 

For the N server queue with all servers having exponential service 
time distributions, but not necessarily identical means, and having a Poisson 
input with parameter X, and having attained equilibrium, the output of each 
individual server will be Poisson if and only if it is equiprobable that a cus- 
tomer be in server i as in server j, for any i,j, conditioned on 1 customers 
in the system. The Poisson outputs will each be independent of all the other 
outputs and will have parameter X 

The proof of this theorem will be presented in two parts, the first part 
showing suffciency and the second part showing necessity. 

For ease of notation and without loss of generality, when referring to 
an individual server, server one with exponential service time parameter 
c j will be chosen. It will be easily seen that extension of what formulas 
are based on referring to an individual server from server one to any other 
server i is only a matter of notational complexity. 


The sufficiency part of the proof follows very closely to Burke's proof. 


( 1 ) 


( 1 ) 


Consequently, only the significant diviations from the work of Burke will 
be given. 
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The necessity part of the proof, will be more involved. The general solu- 
tion of the system of equations will be expanded as a matrix exponential, 
and in order to have a Poisson output various summation equalities will 
result. These equalities will then be shown to necessitate the conditions 
of the theorem. 


V. Sufficiency 

The notation of Burke is redefined to pertain to an individual server 
rather than the general system. Let t^ be any randomly selected instant 
of time, (later taken to be zero for convenience. ) Let L be the length of 
time from t^ to the instant of the next departure from , and define 


F k (t) = Pr [ n(t) = k, L > t - t Q ] 

Note 

oo 

^ F k (t) = F<t) = Pr (L > t - t Q ) 
k=0 

and 


W p k * 

A new symbol, interpreted as the portion of the conditional mean service 
time p contributed by the server j, is defined as 

T| ^ = O' , Q \ 

J J J 


Note 


N 


J = 1 


Before proceeding with the sufficiency part of the proof a lemma is given. 
The lemma is quite comparable to the lemma stated earlier. 


Lemma 2 


If it is equiprobable that a customer be in server i as in server j, for 
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any i, j, conditioned on 1 cuet 9 mers in the system, then the probabilities 
describing the state at any random time selected from the set of departure 
instants of server i are the equilibrium probabilities. 

The conditions of the lemma are simply 


Q- =Qj. v i, j, k 


Thus 


■ 




cr. Q l 

-J 1 


N 




Q j" 


cr. 


The following equations describe the system. 
F 0 (t+dt) =F 0 (t)(l-Xdt) +F 1 (t)(l-Xdt)(|a 1 -i 1 j ) dt 


F.(t+dt) = F._ 1 (t)(Xdt)(l-(i._ jdt) + F.(t)(l-Xdt>(l-|i dt) 

+ F. + 1 (t)(l-Xdt)[ (n l+1 -r, j' 1 ) dt ] 


F N-l (t+dt) =F N-2 (t)(5ldt)(1 - ,,l N-2 dt) +F N . 1 (t)(l-Xdt)(l- l i N _ 1 dt) 
0 ) + F N (t)(l-Xdt)[ (<r - ffl ) dt ] 

F N (t+dt) = F N _ 1 (t)(X dt)(l- n N _ 1 dt) +F N (t)(l-Xdt)(l-<rdt) 

+ F N+1 (t)( 1-X dt) [ < cr - <r 1 ) dt ] 


F (t+dt) = F (t)(Xdt)(l-<rdt) + F.(t)(l-Xdt)(l-<rdt) 
J J “ J 
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+ F j + i (t)( 1- X dt) [ ( a - tr j ) dt ] 

: N +1 < j < CO . 

Performing the usual manipulations, these equations become 

F 0 ' ^ F o + ■ ^ 1 11 1 ^ F 1 


F. = X F. . 

i l-l 


(X+ ^) F i + (^ i+ i- ^ 


i +1 
1 


)F, 


+ 1 


1 < i < N -2 


( 8 ) 


N- 1 


= XF 


N -2 


-(X+u ,)F m . + (or - tr,) F 


N-r n-i 


N 


f' = XF. , - (X+ 0 -) F. + (<r -cr,) F. 7 
J J-l J 1 j+1 


N< j < oo 


This is a system of first order, linear, homogeneous differential equations 

subject to the initial conditions F, (0) = P 7 , where t rt has been chosen as 

k k 0 

zero for convenience. The solution of this system of equations is unique 
and it can be easily verified that with Q* = the solution is 

cr , 

F k (t) =p k e " x ^ * 


The remainder of the proof of the lemma is quite analogous to Burke's 
demonstrations. 

To proceed with the sufficiency part of the theorem's proof some 

symbols are redefined. Let L be the length of time from any departure 

instant of cr^, i. e. , the interdeparture interval length of cr and n(t) be 

the state at time t after some arbitrarily chosen departure instant of n 

t , where for convenience t _ is taken as zero. Define 
d d 
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F k (t) = Pr [ n(t) = k, L > t ] 


and note 
oo 




= F (t) = Pr (L > t) 


With this notation the set of differential equations (8) can be seen to describe 
the system, and with Q . l = Q k the solution to the equations is 


F k (t) = P k e 


-X — t 


Following the format of Burke this is easily used to show that if Qj - Q k , 

v • j 

for all i, j , k, the output is Poisson with parameter A . 

Finally note that by theorem one the overall output is Poisson with 
parameter A and hence has the characteristic function 
- X + X s 

$ (s) = e 

. th 

From the results of this section, the characteristic function of the i out- 

put is ^ 

. - X — 
cr 


f . - X 

(s) = e 


cr, 

+ X — s 

cr 


With 


N 

II ( s 

i=l 1 


N N 

u* cr* 


U * v * 

- x J] — + x Y — 

= e U <s lj cr 


i=l 


i=l 


- X + X s 
= e 

= as \ 

~ \ — / 

it is easy to see from the convolution theorem that all of the outputs are 
independent of each other. 


Necessity 

To prove the necessity part of the theorem, it will be assumed that the 
output from the server is Poisson with parameter [T ♦ 
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The differential equations describing the system may be written in 
matrix notation as 
(10) F' = AF 


with 


A = 


X 

X 

0 


( h -n.) 

-(2+fij) (m 2 -n j 2 ) 

X - ( X + 


0 0 


0 0 


(m 3 -n p 


' U+ V ( ^ k + r n f +1 > 


X -( X +cr) ( r- r j) 


The system of differential equations is subject to the initial conditions 


F (0) = P 

where t^ is a randomly selected instant of time conveniently taken as zero 
and P is the equilibrium probability vector 



It is well known that the solution to these equations may be written as 



= ( I + At + A 2 + * * * ). P 

2 ! 

= P + APt + A 2 P t 2 + . . . 

2 ! 

The product AP can be shown to be 
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AP = 


- X — * P n 

Hi 0 


f+1 

_ x — P 

** 1+1 


- X — ■ P. T , 
cr N- 1 


-\—P. 

cr k 


Since 


£ F (t) = Fr( L > t) 


and the output is Pois son, implying 
Pr(L>t) =e '^1* 

=i - p,t + p^ t z 


2 ! 


Therefore, 


£ = Pi 


and thus it can be shown that 

N-2 


k+1 


— ) P^ 

or k 


p , = x — + x V < — l — 

1 - a r=0 ^k+l 

As before, the characteristic function of the overall output is 


$ (s) = e 


- X + X s 


and the characteristic function of the output of each individual server is 
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Thus again the product is found to be 

n (s) = e" * + X S = $ (s) 
i = l 1 

and therefore the individual outputs must be independent of all the other 
outputs. 

Having established the above result, the remaining portion of the proof 
will proceed under the assumption that the individual outputs are each 
Poisson and independent of the other outputs. 

Combining previous results with the independence assumption gives 

N j. -p j.2 2 

II P [ L> . > t j n{t) = k] -e = 1 - Xt + X + O ( t + ) 

i = 1 

3 

where by O (t +) is meant terms involving the third or higher powers of t. 
Now 

[ Pr L. > t | n (t) = k 1 = p [ L i > n(t) = k ] 

Pr Jn(t) = k] 

and it can be shown that 



84 


A 2 F = 


1 2 


M 1 M ] 

L p 

^ 1^2 

0 

£ +1 

£+2 

^1 

^1 

^£+1 

^£ + 2 


■(A 

\ “i 


2 


P + X' 
e 


£ +1 £+ 2 ' 

’Hi ^1 \ . 

— - i — 1 P + Xp, . 

£ +1 ^ +2 ' 6 * 


£+1 


£ 

*1 




“£+1 ^£ 


' 2 

^ ~T P N-1 + X ^N-ll cr 
cr 

2 

2 ^ 1 

x ~~r p n 

cr 

: 2 

2 ' 4 

x 4- p k 


ff i n i 


N-r 


N-l 


N-l 


AP and A 2 P can be used to form a matrix exponential expansion in 
terms of 0(t 3 +). Continuing this expansion with the previous exponential 
expansion gives 


N k+1 k+2 

y?i rtj 

r . hfi+l ^k+2 

1= 1 


N-l N 

+ 2 E Z 

i- 1 j= i+ 1 


k+1 k+1 

t li 
2 

^k+1 


= 0 < K < N - 3 


and 


N N-l 
r). cr. 

1 1 i_ 

L ll _ _ . cr 


i= 1 


N-l 


N-l IN 

+ 2 E Z 


AN - A IN - 1 

r] . rj : 

LJ = i 


i= 1 j = i + 1 TN-1) 


Note 


( y 4 +1 \ 2 . a fliiV + 2 Y V 'T .- I J 1 
\l?l ^k+J il j l\^k+l/ i=lj=i+l (H k+1) 


1 = 
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Taking the difference gives 


N 


.1 


k+1 k+2 

71 i 2 i 
11 k+1 ^ k+2 




and 


N 


2 

— 


1 




N- 1 
_i 

N-l 


O'. 

1 

O' 



These expressions are most conveniently written as vectors. Let 



Then the expressions become 


X. * 
k 

X k + 1 

X. * X 
k 

X N-1 

* Y s 

X N-I ■* 


k 

X 


N-l 


1 < k < N-2 


where the * indicates the dot product (inner product) of the vectors. 
The vectors are also subject to the constraints 


N 

, 5 , 


k, i 


= 1 


X, . > 0 

k, l - 


N 


y y. = i 

i = 1 1 


Y. > 0 
i — 


The constraints imply the vectors must all terminate on the n 
dimensional hyper plane P defined by the N points 


th 
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P i = ( 8 li * 6 2i ’ * * * * S Ni^ 1 5- 1 1 N 


where is the Kronecker delta. The vector relations may be 
restated as 


x k * < x k + i - x k> = 0 


1 < k < N-2 


N-l * ( Y - X N1 ) = 0 


Thus, either X^ + ^ is same as X^. or the difference X^+j - X^ 
must be perpendicular to X^.» for 1 < k < N-2, and similarly, either 
X ^ ^ is the same as Y or Y-X^ ^ is perpendicular to X^ j . How- 
ever, since each vector terminates on the plane P, their difference 
can be seen as a free vector lying in the plane P. Hence , either 
X^ + i must be the same as X^» or X^ must be perpendicular to 
the plane P . There is only one vector perpendicular to the plane P 
and that is the vector 

x = ( 1/N , 1/N, ... , 1/N) 

Therefore, either X ^ - ^k+1 or ^k = 1 * ^^ 1US » by the above re- 

lations the vectors can be split into two classes. 


^ i • 0 < i < N-l & • X = ± 


X k= Y 


1 < k < i 
i < k < N-l 


Assume i =jt 0. Then 


X = X_ = ... = X. j. 


X i+1 " X i+2 " " X N-1 ” Y 



87 


The following relation can also be shown. 


X 1 P 0 


X 2 P I 


X n-l P N-2 


' X 2 P 0 


X 3 P 1 


X n-1 P N 


+ Y P 


N-2 


For ease of reading the derivation of this relation will not be given until 
after the proof is completed. The above relation may now be rewritten as 


i P n + ..« + j_P. + YP. + YP 

x 0 x i 


i+ 1 


N-2 


P n + . . . + X P. , + YP. + ... YP. 

0 i- 1 i 


N-2 


or 

,R = YP. 

1 x 

or 

x = Y 

Therefore, unless Y = { 1 /N , 1/N, ... , l/N), there is a contradiction 

and i = 0. Hence 


or 


But 


X. = Y Y k 

k 




cr. Q . 

l l 
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Therefore 


or 




Vi, k 


Since the right hand side does not involve the index i, we have Q k = Q k 

i j 

for all i, j, k, thus completing the proof. 

It remains to demonstrate the relation 


X 1 P 0 


X P 
2 1 


X N-l P N-2 


= X 2 P 0 


X P 
3 1 


X N-1 P N - 3 + YP N-2 


The demonstration involves combining the criteria for the outputs to be 
Poisson and the criteria for the outputs to be independent. First the 
criteria for the outputs to be Poisson will be derived. 

As before, define 


F k (t) = Pr [ Lj > t, n(t) = k] 

The system is then described by the system of differential equations 
- F '= A i- F 

subject to 

F(0) = P 

where is the same as earlier 
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The solution is then known to be 

2 

F(t)= e A l* P = (I + Aj ‘ + Aj ^ + .,.)P 

2 t ^ 

= P + AjPt + A i?|l + ... 

2 

where A^P and P are the same as before* 

Now 

^ F k (t) = Pr(L A > t) 
and if the output is Poisson, 

Pr(L 1>t ) = e'Pl^ 1 - Pl t + pJ £r- . 

Thus 

Pr (Lj > t) = V ( P 4 A x Pt + Aj P + . . 

2 

-l-*Pjt+(3i 2 J * • • • 

These results imply 

(J AjP) 2 = j A*P 

Taking the respective sums and comparing leads to 



This is the criteria for the outputs to be Poisson. 
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Now the criteria for the -outputs to be independent will be derived. 
Define 


G^(t) = Pr £ L j ^ 1 1 L»2 ^tj ... » j ^ 1 » n (t) — k ] 


and note 


2^w 


G k (0) = P k 

Pr { L j > tj ^ 1 1 • • • » ^N -> 1 ^ 


The above function is different from the other functions defined earlier 
in that it deals with N-l of the N servers instead of an individual server 
or all the servers. The following equations describe the system. 


° 0 ' " ~ kG 0 + G 1 


i+1 


g . = \G i t - ( \ + M G, + n G. 


V ~i . 11 N l+l 


1 < i < N-2 


(ID 


G N-1 = XG N-2 {x + ^N-l )G N-1 + °N G N 


G' = \G. , 
J # 3-1 


— (X. + <r) G. + o* T G. , 
1 J N 3-I 


N < j < °o 


These equations are subject to the initial conditions 
G (0) = P 

The system of differential equations (11) may be written in matrix 
notation as 


(12) 


G = BG G (0) = P 
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where B is the matrix 



As before, the solution is known to be 

G = e Bt P = (I + Bt + bV + ... )P 

2 1 

= P + BP t + B 2 PT 2 + . . . 

2 ! 

2 

From this result the products BP and B P may be calculated. 
In order to have independence it is necessary that 


Pr(L, > t, > t, ... , , > t) = Pr (L, > t) . , . Pr(L M , >t) 


Thus by associating the t terms of the joint density with the product 


of the densities it is seen to be necessary that 


N- 1 oo 


N-2 oo 


N-l co 


S B 2 P= m Afp ) + 2 l ( l A.P) [l il AP, 

= 0 i= 1 k=0 i=l k=0 *■ j=i+ J. k^O J 
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This leads tQ the requirement 



= 0 

or with sufficient algebra and C 2 may be combined to yield 



and 
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substituted in (13) and with a little manipulation it can be seen that 



But the term on the left is simply for server N implying 



Since this must hold for any server, it follows 
N-2 N-3 

Tj X k+l P k = E X k+2 P k + YP N-2 
k = 0 k= 0 

or 

X 1 P 0 + X 2 P 1 + +X n-l P N-2 


= X 2 P 0 


X 3 P 1 


Y 

N-l N-3 


+ YP 


N-2 


as desired, 
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VIII. Restrictions 

It has been shown that only when G* ~ Q* can the multiple servers 

J K 

with exponential service time distributions having different means each 
generate a Poisson output. It has yet to be shown how this set of equilibrium 
probabilities might be established. In this section the two server case will 
be examined in order to show a means of establishing the equilibrium proba- 
bilities and also the nature of restrictions on the dispersion of the means 
of the service times resulting from the necessity of establishing the equi- 
librium probabilities. 

Consider the two server queue shown in figure two. All the appropriate 
assumptions are made in order to have Poisson outputs from each server. 
Then 



and since their sum must equal one, 

oj = 1 , i = 1. 2 

The following equation describes the system relative to server Qne> wiiere 
R is the probability of assigning a customer arriving to find an empty 
system to server one. 

(14) ( \ + c r_ IP.O 1 = \P.R + a P 

l' l~i Q ' 2~ 2 

The application of the equilibrium analysis given earlier in this paper 



1 - i 

( 1 +^) - l 


shows 



95 


where 


h- 


a l Q l + r 2 Q 2 


P = r p 

1 *1 0 


P = L l P 
2 * *1 0 


Using these results in (14) yields 


(15) 


R = 


o * 1 0- + \ (o-j - or 2 ) 


Thus by assigning customers which arrive to find an empty queue to 
server one with probability R given in (15) each server in the two 
server queue will have a Poisson output. 

Since R is a probability, it is subject to the restriction 


0 < R < 1 


Therefore. 


d cr + \ (o' - o' 9 ) 

0 < —I ^ — < 1 : 

O’ ^ 

c 

Let a j = p <r, then <r 2 = (1 - p) o’ and the following restrictions 

result. ’ ; 


(1 : 6> . i -frir ; ,, j. ; . , 

This restriction is indicated in the table below and is shown in graph 
form in figure three. The^restriction (16) shows that in order to 
establish Poisson outputs from each server the dispersion of the 
service time means must be restricted. 
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£ 

L». B. 


U. B. 

0 

0 


1 

.2 

. 143 


. 781 

.4 

. 222 


. 777 

. 6 

. 273 


. 729 

. 8 

. 308 


. 693 

1. 0 

. 333 


.667 
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FIGURE ONE 


General N Server Queue 
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FIGURE TWO 


Two Server Queue 
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IX. Conclusions 

In this paper it has been shown that a multiple server queue, where 

each server has an exponential service time distribution with possibly 

different means, subject to a Poisson input with parameter X f and after 

attaining equilibrium, will have a Poisson output with parameter X . 

Furthermore, the output of each individual server i will be Ipoisson with 
^i 

parameter X — if and only if it is equiprobable that a customer be in 
server i, conditioned on 1 customers in the system, as in server j, 
with the same conditioning. The two server queue was analyzed to show 
how such equilibrium probabilities might be established and that for such 
equilibrium probabilities to be possible the mean service times must not 
be too different. Since systems are designed by criteria other than mak- 
ing outputs fall into some nice statistical characterization, further in- 
vestigation of the restrictions was not pursued. 
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